	/*This program computes logit demand and simulates corresponding price changes assuming Bertrand competition with product differentiation*/

	#delimit;
	set matsize 10000;
	clear;
	clear matrix;
	mata: mata clear;
	local nbrands=9;
	set seed 070910432;

	*quietly{;
		set mem 250m;
		set more off;
		do "c:\data\restat_12_9\programs\csolve_do\csolve";
		do "c:\data\restat_12_9\programs\csolve_do\callf";
		do "c:\data\restat_12_9\programs\csolve_do\callf_setup";
		do "c:\data\restat_12_9\programs\csolve_do\antilogit_7_8";
		do "c:\data\restat_12_9\programs\csolve_do\anti_costs";
		do "c:\data\restat_12_9\programs\csolve_do\mcosts_antisyr";
		
		
		/*get top level elasticity*/
		mata: mata matuse "c:/data/restat_12_9/analysis_data/top_level";
		
		mata: st_matrix("ols_top_b",ols_b_top[1,1]);
		mata: st_matrix("iv_top_b",iv_b_top[1,1]);
		
		mata: st_matrix("ols_top_mean",ols_b_top);
		mata: st_matrix("iv_top_mean",iv_b_top);
		mata: st_matrix("ols_top_vc",ols_vc_top);
		mata: st_matrix("iv_top_vc",iv_vc_top);

		
		
		scalar ols_top_b=ols_top_b[1,1];
		mata: ols_top_b=ols_b_top[1,1];
		
		scalar iv_top_b=iv_top_b[1,1];
		mata: iv_top_b=iv_b_top[1,1];
		
		*scalar top_b=-1.15;
		
		
		use "c:\data\restat_12_9\analysis_data\hl_syruplogitdata";
		qui tab BRAND2, gen(bra);
		qui tab REGION, gen(reg);
		qui tab t, gen(t);
		qui tab MONTH, gen(m);
		
		local reglist "reg2";
		forval x=3/47{;
			local reglist "`reglist' reg`x'";
		};

		su AVOL;
		scalar vt=r(N)*r(mean);
		forval x=1/`nbrands'{;
			su AVOL if bra`x'==1;
			scalar vshare`x'=r(N)*r(mean)/vt;
		};
		
		
		forval x=1/`nbrands'{;
			su price if bra`x'==1;
			scalar price`x'=r(mean);
		};
		scalar wavgp=vshare1*price1+vshare2*price2+vshare3*price3+vshare4*price4+vshare5*price5+vshare6*price6+vshare7*price7+vshare8*price8+vshare9*price9;
		forval x=1/`nbrands'{;
			forval y=2/47{;
				gen b`x'r`y'=reg`y'*bra`x';
			};
		};
		
		forval x=1/`nbrands'{;
			forval y=2/9{;
				gen m`x'r`y'=bra`x'*m`y';
			};
		};
		
		forval x=1/4{;
			forval y=2/47{;
				local brandxreg "`brandxreg' b`x'r`y'";
			};
		};
		
		forval x=6/7{;
			forval y=2/47{;
				local brandxreg "`brandxreg' b`x'r`y'";
			};
		};
		
		forval x=1/4{;
			forval y=2/9{;
				local brandxmonth "`brandxmonth' m`x'r`y'";
			};
		};
		
		forval x=6/7{;
			forval y=2/9{;
				local brandxmonth "`brandxmonth' m`x'r`y'";
			};
		};
		
		
		
		local logit "price_plprice bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `reglist' t2 t3 t4 t5 t6 t7 t8 t9";
		local logit2 "price_plprice bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `brandxreg' `brandxmonth'";

		
		
		
		
		sort BRAND REGION t;
		by BRAND REGION: gen T=_n;
		sort T REGION BRAND;
		bysort T: gen id=_n;
		tsset id T;
		
		noisily xi: newey2 lshare_plshare `logit'  if BRAND2~="PRIVATE LABEL", lag(1) nocons;
		noisily xi: newey2 lshare_plshare `logit2'  if BRAND2~="PRIVATE LABEL", lag(1) nocons;		
		
		matrix ols_logitB=e(b);
		matrix ols_logitVC=e(V);
		scalar alpha=_b[price_plprice];
		matrix olsb=_b[price_plprice];
		
		
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				sca olse`x'b`y'=-alpha*price`x'*vshare`x'+ols_top_b*vshare`x'*price`x'/wavgp;
			};
		};
		
		forval y=1/`nbrands'{;
			sca olse`y'b`y'=-price`y'*alpha*(1-vshare`y')-ols_top_b*vshare`y'*price`y'/wavgp;
		};
		
		mat olselas=J(9,9,0);
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				mat olselas[`x',`y']=olse`x'b`y';
			};
		};
		mat li olselas;
		
		
		noisily xi: newey2 lshare_plshare bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `reglist' t2 t3 t4 t5 t6 t7 t8 t9 (price_plprice=tiprice1-tiprice9) if BRAND2~="PRIVATE LABEL", lag(1) nocons;
		noisily xi: newey2 lshare_plshare bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `brandxreg' `brandxmonth' (price_plprice=iprice) if BRAND2~="PRIVATE LABEL", lag(1) nocons;
		
		matrix iv_logitB=e(b);
		matrix iv_logitVC=e(V);
		
		scalar ivalpha=_b[price_plprice];
		matrix ivb=_b[price_plprice];
		 
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				sca ive`x'b`y'=-ivalpha*vshare`x'*price`x'+iv_top_b*vshare`x'*price`x'/wavgp;
			};
		};
		forval y=1/`nbrands'{;
			sca ive`y'b`y'=-ivalpha*price`y'*(1-vshare`y')-iv_top_b*vshare`y'*price`y'/wavgp;
		};
		
		mat ivelas=J(9,9,0);
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				mat ivelas[`x',`y']=ive`x'b`y';
			};
		};
		mat li ivelas;
		

		#delimit;

		mat vsbar=J(9,47,0);
		mat pbar=J(9,47,0);
		forval r=1/47{;
			sum AVOL if reg`r'==1;
			sca vt`r'=r(N)*r(mean);
		};
	
	
	
		forval r=1/47{;
			forval b=1/9{;
				qui su AVOL if reg`r'==1&bra`b'==1;
				mat vsbar[`b',`r']=r(N)*r(mean)/vt`r';
				qui su price if reg`r'==1&bra`b'==1;
				mat pbar[`b',`r']=r(mean);
			};	
		};		

		
		
	#delimit cr;
	/*load in post avg prices, and put pre prices and shares in mata.  construct ownership matrix*/
		mata: 
			mata matuse c:/data/restat_12_9/analysis_data/syruphlpost_prices, replace
			
			pbar=st_matrix("pbar")
			vsbar=st_matrix("vsbar")
			postdirect=J(9,47,0)
			postdirect[1,.]=pbar[1,.]*(1+0.003)	/*a.j.*/
			postdirect[2,.]=pbar[2,.]*(1+0.0094)	/*a.j.l.*/
			postdirect[3,.]=pbar[3,.]*(1+0.0160)	/*h.j.*/
			postdirect[4,.]=pbar[4,.]*(1+0.0172)	/*h.j.l.*/
			postdirect[5,.]=pbar[5,.]*(1+0.0323)	/*l.c.*/
			postdirect[6,.]=pbar[6,.]*(1+0.0346)	/*l.c.l.*/
			postdirect[7,.]=pbar[7,.]*(1-0.0090)	/*m.b.*/
			postdirect[8,.]=pbar[8,.]*(1+0.0037)	/*msb light*/
			postdirect[9,.]=pbar[9,.]*(1+0.0072)	/*private label*/
			
			pre_m=I(`nbrands')
			pre_m[1,2]=1
			pre_m[2,1]=1
			pre_m[3,4]=1
			pre_m[4,3]=1
			pre_m[5,6]=1
			pre_m[6,5]=1
			pre_m[7,8]=1
			pre_m[8,7]=1
			pre_m[9,9]=1
			
			
			post_m=pre_m
			post_m[5,7]=1
			post_m[5,8]=1
			post_m[6,7]=1
			post_m[6,8]=1
			
			post_m[7,5]=1
			post_m[7,6]=1
			post_m[8,5]=1
			post_m[8,6]=1
			
			olslogit=J(`nbrands',47,0)
			ivlogit=J(`nbrands',47,0)
			ols_solver_flag=J(47,1,0)
			iv_solver_flag=J(47,1,0)
			olsb=st_matrix("olsb")
			ivb=st_matrix("ivb")
			olsb=-olsb
			ivb=-ivb
			ols_e=-ols_top_b
			iv_e=-iv_top_b
		
			
			for(i=1;i<=cols(pbar);i++){
					pinit=pbar[.,i]
					a=csolve(&anti(),pinit,1e-6,400,vsbar[.,i],pbar[.,i],olsb,ols_e,pre_m,post_m)
					sol=a[1..9,.]
					olslogit[.,i]=sol
					ols_solver_flag[i,1]=a[10,1]
			}
		olslogit=olslogit'
		ans=(olslogit-pbar'):/pbar'
		/*take only markets for which we could compute a solution*/
		ols_sim=100*47/(colsum(J(47,1,1)-!rowsum(ans)))*mean(ans,1)
		
		
		/*Evaluate FOC's*/
		olsfocs=J(47,9,0)
		for(i=1;i<=cols(pbar);i++){
				olsfocs[i,.]=anti(postpbar[.,i],vsbar[.,i],pbar[.,i],olsb,ols_e,pre_m,post_m)'
		}		
		/*couldn't compute solution in some markets*/
		for(i=1;i<=cols(pbar);i++){
					pinit=pbar[.,i]
					a=csolve(&anti(),pinit,1e-6,400,vsbar[.,i],pbar[.,i],ivb,iv_e,pre_m,post_m)
					sol=a[1..9,.]
					ivlogit[.,i]=sol
					iv_solver_flag[i,1]=a[10,1]
		}
		ivlogit=ivlogit'
		ans=(ivlogit-pbar'):/pbar'
		iv_sim=100*47/(colsum(J(47,1,1)-!rowsum(ans)))*mean(ans,1)
		ivfocs=J(47,9,0)
		for(i=1;i<=cols(pbar);i++){
				ivfocs[i,.]=anti(postpbar[.,i],vsbar[.,i],pbar[.,i],ivb,iv_e,pre_m,post_m)'
		}		
		/*couldn't compute solution in some markets*/
		
		
		/*Pre-merger marginal costs, OLS*/
		pre_costs_ols=J(47,9,0)
		for(i=1;i<=cols(pbar);i++){
			pre_costs_ols[i,.]=anti_costs(vsbar[.,i],pbar[.,i],olsb,ols_e,pre_m)'
		}
		
			/*necessary marginal cost changes for observed price increase and simulations to match*/
			post_costs_ols=J(47,9,0)
			for(i=1;i<=cols(pbar);i++){
					init=pre_costs_ols[i,.]'
					a=csolve(&mcosts_antisyr(),init,1e-6,400,vsbar[.,i],postdirect[.,i],olsb,ols_e,pre_m,post_m)
					post_costs_ols[i,.]=a[1..9,.]'				
			}
			mcchange_ols=(post_costs_ols-pre_costs_ols):/pre_costs_ols
			permcchange_ols=100*47/(colsum(J(47,1,1)-!rowsum(mcchange_ols)))*mean(mcchange_ols,1)
		
		/*Pre-merger marginal costs, IV*/
		pre_costs_iv=J(47,9,0)
		for(i=1;i<=cols(pbar);i++){
			pre_costs_iv[i,.]=anti_costs(vsbar[.,i],pbar[.,i],ivb,iv_e,pre_m)'
		}
		
			post_costs_iv=J(47,9,0)
			for(i=1;i<=cols(pbar);i++){
					init=pre_costs_iv[i,.]'
					a=csolve(&mcosts_antisyr(),init,1e-6,400,vsbar[.,i],postdirect[.,i],ivb,iv_e,pre_m,post_m)
					post_costs_iv[i,.]=a[1..9,.]'				
			}
			mcchange_iv=(post_costs_iv-pre_costs_iv):/pre_costs_iv
			permcchange_iv=100*47/(colsum(J(47,1,1)-!rowsum(mcchange_iv)))*mean(mcchange_iv,1)
	end		


/***Now that we have our point estimates of costs and approximate simulated price changes, we draw from asymptotic approx of distributions
of estimators and recalculate for each draw.  will look at empirical distribution of these calculated costs and changes to do inference.****/	
#delimit;
drop _all;


/*number of bootstrap iterations*/
local its=1000;

drawnorm `logit2', n(`its') mean(ols_logitB) cov(ols_logitVC);
mkmat _all;

drop _all;
drawnorm e_topbs b2-b57, n(`its') mean(ols_top_mean) cov(ols_top_vc);
mkmat e_topbs;
drop _all;




use "c:\data\restat_12_9\analysis_data\hl_syruplogitdata";


tab REGION, gen(reg);
tab BRAND2, gen(bra);
forval x=1/47{;
    mat bpbar`x'=J(`nbrands',`its',0);  
	mat bv`x'=J(`nbrands',`its',0);  
    mat bvsbar`x'=J(`nbrands',`its',0);  
};


forval x=1/47{;
    forval w=1/`its'{;
        preserve;
            keep if reg`x'==1;
            bsample;
            forval y=1/`nbrands'{;
                qui sum price if bra`y'==1;
                matrix bpbar`x'[`y',`w']=_result(3);
                qui sum AVOL if bra`y'==1;
                matrix bv`x'[`y',`w']=_result(3);				
            };
        restore;
    };
};
#delimit cr;
mata:
    bp1=st_matrix("bpbar1")
    bp2=st_matrix("bpbar2")
    bp3=st_matrix("bpbar3")
    bp4=st_matrix("bpbar4")
    bp5=st_matrix("bpbar5")
    bp6=st_matrix("bpbar6")
    bp7=st_matrix("bpbar7")
    bp8=st_matrix("bpbar8")
    bp9=st_matrix("bpbar9")
    bp10=st_matrix("bpbar10")
    bp11=st_matrix("bpbar11")
    bp12=st_matrix("bpbar12")
    bp13=st_matrix("bpbar13")
    bp14=st_matrix("bpbar14")
    bp15=st_matrix("bpbar15")
    bp16=st_matrix("bpbar16")
    bp17=st_matrix("bpbar17")
    bp18=st_matrix("bpbar18")
    bp19=st_matrix("bpbar19")
    bp20=st_matrix("bpbar20")
    bp21=st_matrix("bpbar21")
    bp22=st_matrix("bpbar22")
    bp23=st_matrix("bpbar23")
    bp24=st_matrix("bpbar24")
    bp25=st_matrix("bpbar25")
    bp26=st_matrix("bpbar26")
    bp27=st_matrix("bpbar27")
    bp28=st_matrix("bpbar28")
    bp29=st_matrix("bpbar29")
    bp30=st_matrix("bpbar30")
    bp31=st_matrix("bpbar31")
    bp32=st_matrix("bpbar32")
    bp33=st_matrix("bpbar33")
    bp34=st_matrix("bpbar34")
    bp35=st_matrix("bpbar35")
    bp36=st_matrix("bpbar36")
    bp37=st_matrix("bpbar37")
    bp38=st_matrix("bpbar38")
    bp39=st_matrix("bpbar39")
    bp40=st_matrix("bpbar40")
    bp41=st_matrix("bpbar41")
    bp42=st_matrix("bpbar42")
    bp43=st_matrix("bpbar43")
    bp44=st_matrix("bpbar44")
    bp45=st_matrix("bpbar45")
    bp46=st_matrix("bpbar46")
    bp47=st_matrix("bpbar47")

    ppoint=J(1,47,NULL)
    ppoint[1]=&bp1
    ppoint[2]=&bp2
    ppoint[3]=&bp3
    ppoint[4]=&bp4
    ppoint[5]=&bp5
    ppoint[6]=&bp6
    ppoint[7]=&bp7
    ppoint[8]=&bp8
    ppoint[9]=&bp9
    ppoint[10]=&bp10
    ppoint[11]=&bp11
    ppoint[12]=&bp12
    ppoint[13]=&bp13
    ppoint[14]=&bp14
    ppoint[15]=&bp15
    ppoint[16]=&bp16
    ppoint[17]=&bp17
    ppoint[18]=&bp18
    ppoint[19]=&bp19
    ppoint[20]=&bp20
    ppoint[21]=&bp21
    ppoint[22]=&bp22
    ppoint[23]=&bp23
    ppoint[24]=&bp24
    ppoint[25]=&bp25
    ppoint[26]=&bp26
    ppoint[27]=&bp27
    ppoint[28]=&bp28
    ppoint[29]=&bp29
    ppoint[30]=&bp30
    ppoint[31]=&bp31
    ppoint[32]=&bp32
    ppoint[33]=&bp33
    ppoint[34]=&bp34
    ppoint[35]=&bp35
    ppoint[36]=&bp36
    ppoint[37]=&bp37
    ppoint[38]=&bp38
    ppoint[39]=&bp39
    ppoint[40]=&bp40
    ppoint[41]=&bp41
    ppoint[42]=&bp42
    ppoint[43]=&bp43
    ppoint[44]=&bp44
    ppoint[45]=&bp45
    ppoint[46]=&bp46
    ppoint[47]=&bp47    

    bv1=st_matrix("bv1")
    bv2=st_matrix("bv2")
    bv3=st_matrix("bv3")
    bv4=st_matrix("bv4")
    bv5=st_matrix("bv5")
    bv6=st_matrix("bv6")
    bv7=st_matrix("bv7")
    bv8=st_matrix("bv8")
    bv9=st_matrix("bv9")
    bv10=st_matrix("bv10")
    bv11=st_matrix("bv11")
    bv12=st_matrix("bv12")
    bv13=st_matrix("bv13")
    bv14=st_matrix("bv14")
    bv15=st_matrix("bv15")
    bv16=st_matrix("bv16")
    bv17=st_matrix("bv17")
    bv18=st_matrix("bv18")
    bv19=st_matrix("bv19")
    bv20=st_matrix("bv20")
    bv21=st_matrix("bv21")
    bv22=st_matrix("bv22")
    bv23=st_matrix("bv23")
    bv24=st_matrix("bv24")
    bv25=st_matrix("bv25")
    bv26=st_matrix("bv26")
    bv27=st_matrix("bv27")
    bv28=st_matrix("bv28")
    bv29=st_matrix("bv29")
    bv30=st_matrix("bv30")
    bv31=st_matrix("bv31")
    bv32=st_matrix("bv32")
    bv33=st_matrix("bv33")
    bv34=st_matrix("bv34")
    bv35=st_matrix("bv35")
    bv36=st_matrix("bv36")
    bv37=st_matrix("bv37")
    bv38=st_matrix("bv38")
    bv39=st_matrix("bv39")
    bv40=st_matrix("bv40")
    bv41=st_matrix("bv41")
    bv42=st_matrix("bv42")
    bv43=st_matrix("bv43")
    bv44=st_matrix("bv44")
    bv45=st_matrix("bv45")
    bv46=st_matrix("bv46")
    bv47=st_matrix("bv47")
    vpoint=J(1,47,NULL)
    vpoint[1]=&bv1
    vpoint[2]=&bv2
    vpoint[3]=&bv3
    vpoint[4]=&bv4
    vpoint[5]=&bv5
    vpoint[6]=&bv6
    vpoint[7]=&bv7
    vpoint[8]=&bv8
    vpoint[9]=&bv9
    vpoint[10]=&bv10
    vpoint[11]=&bv11
    vpoint[12]=&bv12
    vpoint[13]=&bv13
    vpoint[14]=&bv14
    vpoint[15]=&bv15
    vpoint[16]=&bv16
    vpoint[17]=&bv17
    vpoint[18]=&bv18
    vpoint[19]=&bv19
    vpoint[20]=&bv20
    vpoint[21]=&bv21
    vpoint[22]=&bv22
    vpoint[23]=&bv23
    vpoint[24]=&bv24
    vpoint[25]=&bv25
    vpoint[26]=&bv26
    vpoint[27]=&bv27
    vpoint[28]=&bv28
    vpoint[29]=&bv29
    vpoint[30]=&bv30
    vpoint[31]=&bv31
    vpoint[32]=&bv32
    vpoint[33]=&bv33
    vpoint[34]=&bv34
    vpoint[35]=&bv35
    vpoint[36]=&bv36
    vpoint[37]=&bv37
    vpoint[38]=&bv38
    vpoint[39]=&bv39
    vpoint[40]=&bv40
    vpoint[41]=&bv41
    vpoint[42]=&bv42
    vpoint[43]=&bv43
    vpoint[44]=&bv44
    vpoint[45]=&bv45
    vpoint[46]=&bv46
    vpoint[47]=&bv47

    b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")

end
    
	
	
	
set more off	
/*OLS Elasticities Here*/
mata:
    e1j=J(1000,9,0)
    e2j=J(1000,9,0)
    e3j=J(1000,9,0)
    e4j=J(1000,9,0)
    e5j=J(1000,9,0)
    e6j=J(1000,9,0)
    e7j=J(1000,9,0)
    e8j=J(1000,9,0)
    e9j=J(1000,9,0)
	

for(z=1;z<=1000;z++){
        /*makes price coefficient matrix for each draw*/
	bz=-b[z,1]
	ez=-e[z,1]
	

	/*need to get vs, p*/
	vx=J(9,47,0)
	px=J(9,47,0)

	for(j=1;j<=47;j++){           
		px[.,j]=(*ppoint[j])[.,z]
		vx[.,j]=(*vpoint[j])[.,z]
	}	            

	v=(rowsum(vx):/colsum(rowsum(vx)))'
	p=mean(px',1)
	
	elas=J(9,9,0)
	wavgp=v*p'
	for(m=1;m<=9;m++){
		for(n=1;n<=9;n++){
			elas[m,n]=bz*v[m]*p[m]-ez*v[m]*p[m]/wavgp
		}
	}
	
	for(m=1;m<=9;m++){
		elas[m,m]=bz*p[m]*(1-v[m])-ez*v[m]*p[m]/wavgp
	}
		



e1j[z,1]=elas[1,1]
e1j[z,2]=elas[1,2]
e1j[z,3]=elas[1,3]
e1j[z,4]=elas[1,4]
e1j[z,5]=elas[1,5]
e1j[z,6]=elas[1,6]
e1j[z,7]=elas[1,7]
e1j[z,8]=elas[1,8]
e1j[z,9]=elas[1,9]

e2j[z,1]=elas[2,1]
e2j[z,2]=elas[2,2]
e2j[z,3]=elas[2,3]
e2j[z,4]=elas[2,4]
e2j[z,5]=elas[2,5]
e2j[z,6]=elas[2,6]
e2j[z,7]=elas[2,7]
e2j[z,8]=elas[2,8]
e2j[z,9]=elas[2,9]

e3j[z,1]=elas[3,1]
e3j[z,2]=elas[3,2]
e3j[z,3]=elas[3,3]
e3j[z,4]=elas[3,4]
e3j[z,5]=elas[3,5]
e3j[z,6]=elas[3,6]
e3j[z,7]=elas[3,7]
e3j[z,8]=elas[3,8]
e3j[z,9]=elas[3,9]

e4j[z,1]=elas[4,1]
e4j[z,2]=elas[4,2]
e4j[z,3]=elas[4,3]
e4j[z,4]=elas[4,4]
e4j[z,5]=elas[4,5]
e4j[z,6]=elas[4,6]
e4j[z,7]=elas[4,7]
e4j[z,8]=elas[4,8]
e4j[z,9]=elas[4,9]

e5j[z,1]=elas[5,1]
e5j[z,2]=elas[5,2]
e5j[z,3]=elas[5,3]
e5j[z,4]=elas[5,4]
e5j[z,5]=elas[5,5]
e5j[z,6]=elas[5,6]
e5j[z,7]=elas[5,7]
e5j[z,8]=elas[5,8]
e5j[z,9]=elas[5,9]

e6j[z,1]=elas[6,1]
e6j[z,2]=elas[6,2]
e6j[z,3]=elas[6,3]
e6j[z,4]=elas[6,4]
e6j[z,5]=elas[6,5]
e6j[z,6]=elas[6,6]
e6j[z,7]=elas[6,7]
e6j[z,8]=elas[6,8]
e6j[z,9]=elas[6,9]

e7j[z,1]=elas[7,1]
e7j[z,2]=elas[7,2]
e7j[z,3]=elas[7,3]
e7j[z,4]=elas[7,4]
e7j[z,5]=elas[7,5]
e7j[z,6]=elas[7,6]
e7j[z,7]=elas[7,7]
e7j[z,8]=elas[7,8]
e7j[z,9]=elas[7,9]

e8j[z,1]=elas[8,1]
e8j[z,2]=elas[8,2]
e8j[z,3]=elas[8,3]
e8j[z,4]=elas[8,4]
e8j[z,5]=elas[8,5]
e8j[z,6]=elas[8,6]
e8j[z,7]=elas[8,7]
e8j[z,8]=elas[8,8]
e8j[z,9]=elas[8,9]

e9j[z,1]=elas[9,1]
e9j[z,2]=elas[9,2]
e9j[z,3]=elas[9,3]
e9j[z,4]=elas[9,4]
e9j[z,5]=elas[9,5]
e9j[z,6]=elas[9,6]
e9j[z,7]=elas[9,7]
e9j[z,8]=elas[9,8]
e9j[z,9]=elas[9,9]
}

se_elas=J(9,9,0)

se_elas[1,.]=sqrt(mean(e1j:*e1j)-mean(e1j):*mean(e1j))				      
se_elas[2,.]=sqrt(mean(e2j:*e2j)-mean(e2j):*mean(e2j))				      
se_elas[3,.]=sqrt(mean(e3j:*e3j)-mean(e3j):*mean(e3j))
se_elas[4,.]=sqrt(mean(e4j:*e4j)-mean(e4j):*mean(e4j))
se_elas[5,.]=sqrt(mean(e5j:*e5j)-mean(e5j):*mean(e5j))				      
se_elas[6,.]=sqrt(mean(e6j:*e6j)-mean(e6j):*mean(e6j))				      
se_elas[7,.]=sqrt(mean(e7j:*e7j)-mean(e7j):*mean(e7j))				      
se_elas[8,.]=sqrt(mean(e8j:*e8j)-mean(e8j):*mean(e8j))				      
se_elas[9,.]=sqrt(mean(e9j:*e9j)-mean(e9j):*mean(e9j))				      				      

end	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
local nbrands=9	
mata:	
	bslogit=J(1000,`nbrands',0)
	bslogitfocs=J(1000,`nbrands',0)
	for(z=1;z<=1000;z++){
		bz=-b[z,1]
		ez=-e[z,1]		
        temp=J(47,9,0)

		for(j=1;j<=47;j++){
            	p=(*ppoint[j])[.,z]
	            pinit=(*ppoint[j])[.,z]
				v=(*vpoint[j])[.,z]
				vsbar=v:/colsum(v)
				
	            a=csolve(&anti(),pinit,1e-6,400,vsbar,p,bz,ez,pre_m,post_m)
      	        
				a=a[1..9,.]
            	a=(a-p):/p
	            temp[j,.]=a'
	      }
		
		
		bslogit[z,.]=100*47/(colsum(J(47,1,1)-!rowsum(temp)))*mean(temp,1)
	
		olsfocs=J(47,9,0)
		for(j=1;j<=47;j++){
		    	p=(*ppoint[j])[.,z]	            
	            v=(*vpoint[j])[.,z]
				vsbar=v:/colsum(v)				
				olsfocs[j,.]=anti(postdirect[.,j],vsbar,p,bz,ez,pre_m,post_m)'
		}	
	bslogitfocs[z,.]=mean(olsfocs)
	}
end
mata: st_matrix("bslogit", bslogit)
mata: st_matrix("bslogitfocs", bslogitfocs)
drop _all
svmat bslogit 
svmat bslogitfocs
log using syrtable2col9, replace
centile bslogit1, centile(2.5, 97.5)
centile bslogit2, centile(2.5, 97.5)
centile bslogit3, centile(2.5, 97.5)
centile bslogit4, centile(2.5, 97.5)
centile bslogit5, centile(2.5, 97.5)
centile bslogit6, centile(2.5, 97.5)
centile bslogit7, centile(2.5, 97.5)

centile bslogitfocs1, centile(2.5, 97.5)
centile bslogitfocs2, centile(2.5, 97.5)
centile bslogitfocs3, centile(2.5, 97.5)
centile bslogitfocs4, centile(2.5, 97.5)
centile bslogitfocs5, centile(2.5, 97.5)
centile bslogitfocs6, centile(2.5, 97.5)
centile bslogitfocs7, centile(2.5, 97.5)
centile bslogitfocs8, centile(2.5, 97.5)
centile bslogitfocs9, centile(2.5, 97.5)
log close


save "c:\data\restat_12_9\analysis_data\ols_hlsyruplogit", replace






#delimit;
set more off;
drop _all;
/*number of bootstrap iterations*/

set seed 070910432;

local its=1000;
	local reglist "reg2";
	forval x=3/47{;
		local reglist "`reglist' reg`x'";
	};


	
local logit "price_plprice bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `brandxreg' `brandxmonth'";
drawnorm `logit', n(`its') mean(iv_logitB) cov(iv_logitVC);
mkmat _all;
drop _all;
drawnorm e_topbs b2-b57, n(`its') mean(iv_top_mean) cov(iv_top_vc);
mkmat e_topbs;
drop _all;


#delimit cr
/*IV Elasticities Here*/
mata:
	b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")
    e1j=J(1000,9,0)
    e2j=J(1000,9,0)
    e3j=J(1000,9,0)
    e4j=J(1000,9,0)
    e5j=J(1000,9,0)
    e6j=J(1000,9,0)
    e7j=J(1000,9,0)
    e8j=J(1000,9,0)
    e9j=J(1000,9,0)
	

for(z=1;z<=1000;z++){
        /*makes price coefficient matrix for each draw*/
	bz=-b[z,1]
	ez=-e[z,1]
	
	/*need to get vs, p*/
	vx=J(9,47,0)
	px=J(9,47,0)

	for(j=1;j<=47;j++){           
		px[.,j]=(*ppoint[j])[.,z]
		vx[.,j]=(*vpoint[j])[.,z]
	}	            

	v=(rowsum(vx):/colsum(rowsum(vx)))'
	p=mean(px',1)

	
	elas=J(9,9,0)
	wavgp=v*p'
	for(m=1;m<=9;m++){
		for(n=1;n<=9;n++){
			elas[m,n]=bz*v[m]*p[m]-ez*v[m]*p[m]/wavgp
		}
	}
	
	for(m=1;m<=9;m++){
		elas[m,m]=bz*p[m]*(1-v[m])-ez*v[m]*p[m]/wavgp
	}
		



e1j[z,1]=elas[1,1]
e1j[z,2]=elas[1,2]
e1j[z,3]=elas[1,3]
e1j[z,4]=elas[1,4]
e1j[z,5]=elas[1,5]
e1j[z,6]=elas[1,6]
e1j[z,7]=elas[1,7]
e1j[z,8]=elas[1,8]
e1j[z,9]=elas[1,9]

e2j[z,1]=elas[2,1]
e2j[z,2]=elas[2,2]
e2j[z,3]=elas[2,3]
e2j[z,4]=elas[2,4]
e2j[z,5]=elas[2,5]
e2j[z,6]=elas[2,6]
e2j[z,7]=elas[2,7]
e2j[z,8]=elas[2,8]
e2j[z,9]=elas[2,9]

e3j[z,1]=elas[3,1]
e3j[z,2]=elas[3,2]
e3j[z,3]=elas[3,3]
e3j[z,4]=elas[3,4]
e3j[z,5]=elas[3,5]
e3j[z,6]=elas[3,6]
e3j[z,7]=elas[3,7]
e3j[z,8]=elas[3,8]
e3j[z,9]=elas[3,9]

e4j[z,1]=elas[4,1]
e4j[z,2]=elas[4,2]
e4j[z,3]=elas[4,3]
e4j[z,4]=elas[4,4]
e4j[z,5]=elas[4,5]
e4j[z,6]=elas[4,6]
e4j[z,7]=elas[4,7]
e4j[z,8]=elas[4,8]
e4j[z,9]=elas[4,9]

e5j[z,1]=elas[5,1]
e5j[z,2]=elas[5,2]
e5j[z,3]=elas[5,3]
e5j[z,4]=elas[5,4]
e5j[z,5]=elas[5,5]
e5j[z,6]=elas[5,6]
e5j[z,7]=elas[5,7]
e5j[z,8]=elas[5,8]
e5j[z,9]=elas[5,9]

e6j[z,1]=elas[6,1]
e6j[z,2]=elas[6,2]
e6j[z,3]=elas[6,3]
e6j[z,4]=elas[6,4]
e6j[z,5]=elas[6,5]
e6j[z,6]=elas[6,6]
e6j[z,7]=elas[6,7]
e6j[z,8]=elas[6,8]
e6j[z,9]=elas[6,9]

e7j[z,1]=elas[7,1]
e7j[z,2]=elas[7,2]
e7j[z,3]=elas[7,3]
e7j[z,4]=elas[7,4]
e7j[z,5]=elas[7,5]
e7j[z,6]=elas[7,6]
e7j[z,7]=elas[7,7]
e7j[z,8]=elas[7,8]
e7j[z,9]=elas[7,9]

e8j[z,1]=elas[8,1]
e8j[z,2]=elas[8,2]
e8j[z,3]=elas[8,3]
e8j[z,4]=elas[8,4]
e8j[z,5]=elas[8,5]
e8j[z,6]=elas[8,6]
e8j[z,7]=elas[8,7]
e8j[z,8]=elas[8,8]
e8j[z,9]=elas[8,9]

e9j[z,1]=elas[9,1]
e9j[z,2]=elas[9,2]
e9j[z,3]=elas[9,3]
e9j[z,4]=elas[9,4]
e9j[z,5]=elas[9,5]
e9j[z,6]=elas[9,6]
e9j[z,7]=elas[9,7]
e9j[z,8]=elas[9,8]
e9j[z,9]=elas[9,9]
}

se_elas=J(9,9,0)

se_elas[1,.]=sqrt(mean(e1j:*e1j)-mean(e1j):*mean(e1j))				      
se_elas[2,.]=sqrt(mean(e2j:*e2j)-mean(e2j):*mean(e2j))				      
se_elas[3,.]=sqrt(mean(e3j:*e3j)-mean(e3j):*mean(e3j))
se_elas[4,.]=sqrt(mean(e4j:*e4j)-mean(e4j):*mean(e4j))
se_elas[5,.]=sqrt(mean(e5j:*e5j)-mean(e5j):*mean(e5j))				      
se_elas[6,.]=sqrt(mean(e6j:*e6j)-mean(e6j):*mean(e6j))				      
se_elas[7,.]=sqrt(mean(e7j:*e7j)-mean(e7j):*mean(e7j))				      
se_elas[8,.]=sqrt(mean(e8j:*e8j)-mean(e8j):*mean(e8j))				      
se_elas[9,.]=sqrt(mean(e9j:*e9j)-mean(e9j):*mean(e9j))				      				      

end	




set more off
local nbrands=9
mata:
    b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")
    bslogit=J(1000,9,0)
	bsivfocs=J(1000,`nbrands',0)

	for(z=1;z<=1000;z++){
		bz=-b[z,1]
		ez=-e[z,1]		
        temp=J(47,9,0)

		for(j=1;j<=47;j++){
            	p=(*ppoint[j])[.,z]
	            pinit=(*ppoint[j])[.,z]
				v=(*vpoint[j])[.,z]
				vsbar=v:/colsum(v)
				
	            a=csolve(&anti(),pinit,1e-6,400,vsbar,p,bz,ez,pre_m,post_m)
      	        
				a=a[1..9,.]
            	a=(a-p):/p
	            temp[j,.]=a'
	      }
		
		
		bslogit[z,.]=100*47/(colsum(J(47,1,1)-!rowsum(temp)))*mean(temp,1)
	
		ivfocs=J(47,9,0)
		for(j=1;j<=47;j++){
		    	p=(*ppoint[j])[.,z]	            
	            v=(*vpoint[j])[.,z]
				vsbar=v:/colsum(v)				
				ivfocs[j,.]=anti(postdirect[.,j],vsbar,p,bz,ez,pre_m,post_m)'
		}	
	bsivfocs[z,.]=mean(ivfocs)
	}
end
mata: st_matrix("bslogit", bslogit)
mata: st_matrix("bsivfocs", bsivfocs)
drop _all
svmat bslogit
svmat bsivfocs
log using syrtable2col10, replace
centile bslogit1, centile(2.5, 97.5)
centile bslogit2, centile(2.5, 97.5)
centile bslogit3, centile(2.5, 97.5)
centile bslogit4, centile(2.5, 97.5)
centile bslogit5, centile(2.5, 97.5)
centile bslogit6, centile(2.5, 97.5)
centile bslogit7, centile(2.5, 97.5)

centile bsivfocs1, centile(2.5, 97.5)
centile bsivfocs2, centile(2.5, 97.5)
centile bsivfocs3, centile(2.5, 97.5)
centile bsivfocs4, centile(2.5, 97.5)
centile bsivfocs5, centile(2.5, 97.5)
centile bsivfocs6, centile(2.5, 97.5)
centile bsivfocs7, centile(2.5, 97.5)
centile bsivfocs8, centile(2.5, 97.5)
centile bsivfocs9, centile(2.5, 97.5)

log close


save "c:\data\restat_12_9\analysis_data\iv_hlsyruplogit", replace
